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Table Captions 


Table i: Normalized RMS Model Errors in Streamf unction and Velocity 
fu,v) After One Period for Njnlinear Topographic Wave Solutions. 


Table II: Required cpu Time Per Model Time-Step on a VAXll-780. 


ABSTRACT 

A prototype four-dimensional (x.y.z.t) continental shelf /deep ocean 
model is described. In its present form, the model incorporates the ef- 
fects of finite-amplitude topography, advective nonlinearities, and vari- 
able stratification and rotation. The model can be forced either directly 
by imposed atmospheric windstress and surface pressure distributions, and 
energetic mean currents imposed by the exterior oceanic circulation; or 
indirectly by initial distributions of shoreward propagating mesoscale 
waves and eddies. 

To avoid concerns over the appropriate specification of "open" bound- 
ary conditions on the cross-shelf and seaward model boundaries, a peri- 
odic channel geometry (oriented along-coast) is used. The model employs 
a traditional finite-difference expansion in the cross-shelf direction, 
and a Fourier (periodic) representation in the long-shelf coordinate. 

A modified sigma coordinate system, and a Chebyshev-tau approximation 
scheme, are used to incorporate the vertical dependence. 

The model has been validated against a variety of propagating top- 
ographic wave problems. Representative run times and error estimates 
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1. Introduction 

The dynamical coupling between the large-scale ocean circulation and 
motions on the continental shelf and slope is at present poorly under- 
stood; nonetheless, the existence and implications of that coupling are 
of concern to both coastal and deep-water oceanographers. For example, 
exterior current systems undoubtedly place strong constraints on the form 
of the adjacent shelf circulation, thereby dictating circulation patterns 
in the important near-shore zone. Increasing evidence suggests, for in- 
stance, that the observed southwesterly directed mean flow along the 
northeastern continental shelf of the United States is driven by a long- 
shore pressure gradient supported by the large-scale ocean circulation 
(Csanady, 1978). Conversely, flow along the continental margin can be 
thought of as a boundary layer component of the exterior ocean circula- 
tion (Beardsley and Winant, 1979). A refined description of the- closure 
of open ocean flow by the shelf/slope system would therefore contribute 
strongly to modeling studies of the wind- and eddy-driven oceanic general 
circulation [see also, Haidvogel (1979)J. 

These and other considerations suggest that a detailed knowledge of 
the operative continental shelf/deep ocean coupling mechanisms is of con- 
siderable practical importance. This is not only true, however, from 
strict, y an ocean modeling point of view. As well as being areas of in- 
tense mean flow and potentially significant mesoscale eddy/mean field in- 
teraction, the western boundaries of the world's oceans are also sites of 
strong ocean-atmosphere transfers - for example, of heat. The continental 
slope regions - representing as they do a large reservoir of heat, and 
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acting as preferred locations for oceanic frontal formation - may con- 
tribute strongly to the maintenance and form of that heat flux. A proper 
definition of the role of the continental margins in closing the general 
ocean circulation is ultimately relevant, therefore, to a broad range of 
questions, including local oceanographic and global air-sea interaction 
(i.e., climate) problems. 

Recent observational evidence in the near-shore and deep ocean re- 
gions suggest three plausible mechanisms for shelf/oper ocean coupling. 

The first is the aforementioned mean alongshore pressure gradient (be- 
lieved to be) imposed at the outer edge of the shelf by the exterior 
large-scale circulation. Such a pressure gradient, and its extension 
across the narrow shelf, is thought to drive the mean southwestward flow 
in the Middle Atlantic Bight and thereby to support the observed shelf- 
slope front (Beardsley and Winant, 1979). Second are nearly linear, low 
frequency motions — thought to be topographic Rossby waves — which, gener- 
ated in the deeper ocean, are directed shoreward and propagate onto the 
continental slope. The most complete description of the properties of 
these low frequency motions has been synthesized from data taken at 
Site D, a long-term observing location situated on the continental rise 
north of the mean axis of the Gulf Stream. Results at Site D indicate 
most of the low frequency energy to be associated with periods of a week 
to a few months and with a shoreward energy flux. Thompson (1977) offers 
a review. Finally, mixing and/or transport of properties across the shelf- 
slope boundary is associated with intense eddy motions impinging on the 
shelf/slope region from the open ocean. These nonlinear features, though 


intermittent, are strongly energetic and have been observed to account 
for sizable local heat and momentum exchange [e.g., see Smith {in press)]. 
Such intense mesoscale features can often be seen, and their penetration 
of shelf waters monitored, by satellite infrared imagery (Ring Workshop, 
1977). 

Theoretical and numerical studies of linear topographic Rossby waves 
incident from the open ocean onto the continental slope have been made 
(Kroll and Niiler, 1976; Ou, 1979). Results suggest that linear wave 
mechanisms are unlikely to be quantitatively important in the coupling 
of shelf and open ocean circulations. At the present time, therefore, 
the "best-guess" scenario is that: (1) longshore pressure gradients im- 
posed by the exterior circulation are instrumental in supporting mean 
shelf currents; and (2) energetic eddy motions provide significant 
transient exchange of heat and momentum across the continental slope. 

A numerical model, using primitive equation dynamics, has been devel- 
oped to begin to investigate the details of this continental shelf/deep 
ocean coupling. The dynamic equations and physical assumptions of this 
model are quickly reviewed in Section 2. A sigma coordinate system for 
the vertical dependence is then introduced (Section 3), and the resulting 
four-dimensional problem (in x,y,a and t) is discretized for efficient 
and accurate solution (Section 4). Test problems are examined in Sec- 
tion 5 to determine model accuracy and cpu requirements. 

2. Model Equations 


The model is governed by the hydrostatic primitive equations, which 
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can be written: 
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where, in standard notation: 

(u,v,w) = (x,y,z) components of vector velocity v 


p = density 



imposed (atmospheric) pressure 

V 

[1 e(y)] = Coriolis parameter 

acceleration of gravity (acting vertically downwards) 

vertical diffusive coefficients for horizontal 
momentum and density. 


Equations (1) and (2) express the momentum balance in the x and y direc- 
tions, respectively. The time evolution of the density field, p(x,y,z,t), 
is go/erned by the advecti ve-diffusive equation (3). Lastly, equations (4) 
and (5) express the hydrostaticity and incompressibility of the fr:>w field. 
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Note that the Boussinesq and rigid lid approximations have been invoked. 
The latter, in particular, is numerically advantageous in that it effec- 
tively filters out surface gravity waves which otherwise place strong 
restrictions on time-stepping (see Section 4), We do not expect, however, 
that this assumption will seriously affect the physical representation of 
the coastal/deep ocean coupling of interest. 

For convenience, and to avoid dealing with the considerable problem 
associated with the specification of appropriate "open" boundary condi- 
tions, we choose to work within a periodic channel oriented in the along- 
coast direction. The region is bounded above by a rigid lid, and below 
by a variable bottom of depth h(x,y). The assumption cF periodicity in 
the along-coast direction, and the placement of a fictitious "wall" at 
some distance from the coast are, of course, only crude approximations 
to physical reality. Although they allow rather confident and straight- 
forward treatment of equation (1-5), the influence of these assumptions 
on the modeled flow must be borne in mind. It is expected, however, that 
the existence of these boundaries will be unimportant if they are placed 
at a sufficient distance from the sites of coastal/deep ocean interaction. 

Consistent with the assumption of the periodic channel geometry, suit- 
able boundary conditions for equations (1-5) are chosen as follows: 
sidewalls (y = 0, L^) v = 0 (6) 


top (z = 0) 
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or u = 0 

or V = 0 (8) 

Thus, the sidewalls are required to be impermeable (equation 6), and sur- 
face distributions of wind stress and density flux (t^) are 

prescribed on the top (equation 7). On the variable bottom, z » -h(x,y), 
the horizontal velocity components are constrained either to vanish, or 
to accommodate a prescribed bottom stress; the density flux is taken to 
vanish. It is relevant to note at this point that the vertical (Cheby- 
shev-tau) discretization scheme to be described in next section can triv- 
ially accommodate quite arbitrary boundary conditions on z = 0, -h. 

Lastly, all variables are required to be periodic over the interval 
0 < X < L^. 

3. Sigma (stretched vertical) Coordinate System 

From the point of view of the computational model, it is highly con- 
venient to introduce a stretched vertical coordinate system, which essen- 
tially "flattens out" the variable bottom at z = -h(x,y). Such "sigma" 
coordinate systems have long been used, with slight appropriate modifi- 
cation, in both meteorology and oceanography [e.g., Phillips (1957) and 
Freeman et al. (1972)]. To proceed, we make the coordinate transformation: 

A 

X = x 

y = y 

0 = 1+2 (z/h) 


and bottom (z = -h) v = 


= T ; (x.y) 


3Z 

11 
3Z 

^ = 0 . 
3Z 
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and t = t . 

In the stretched system, the vertical coordinate a spans the range -1 £ o 
< 1; we are therefore left with a domain of regular shape. As a trade-off 
for this geometric simplification, however, the dynamic equations become 
somewhat more complicated. The resulting dynamic equations are, dropping 
the carats: 

- fv = - V . vu - |_ [ (h p do) - (1 - a) p ] 

( 10 ) 
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ifChu] * I^^Chv] * h|S = 0 (13) 

where the vertical velocity in sigma coordinates: 

w (x,y,o,t) = ~ [(l-a)u + (l-o) V |y + 2w(x,y,z,t)] . 

p° is the net pressure at z = 0 - the sum of the imposed atmospheric 
pressure and the "reaction pressure" due to the presence of the rigid lid. 
Note also that the dynamic pressure has been eliminated from the equations 
using the hydrostatic relation. 
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In the stretched coordinate, boundary conditions (6-8) become: 

V « 0 


sidewalls (y = 0, L^) 


top (a = 1) 


(M? ' ’ X (>'•» 
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and bottom (a = -1) 


x 2V n 3U 
^h ^ 3a 

(li) ii 

'h ' 30 

/ 2n \ 3p 
^h ^ 3a 


= T (x.y) 


- T ^ (x,y) 


= 0 . 


or u = 0 
or V = 0 


( 16 ) 


Details of this coordinate transformation have been given by several 
authors - see, for instance, Owen (1980). 


4. Numerical Discretization and Solution 
(4a) Chebyshev Polynomial Basis Set 

In the vertical (a) direction, it is dynamically appropriate and 
computationally convenient to represent the motion in terms of a sum of 
its depth-averaged (external) and depth-varying (internal) components. 

To do so, we represent the vertical dependence of the dependent variables 
as an expansion in the polynomial set P|^(o): 


(u,v,p,w) 


N 

I CU|^(x,y,t), V|^(x,y,t), p,^(x,y,t), W|^(x,y,t)] P,^(o). 
k“0 


(17) 


We further require that 


/^lPl^(o) da = 6 


Ok 


- that is, that only the zeroeth 
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order polynomial have a non-zero vertical integral. 

The Pj^(o) expansion functions can be chosen somewhat arbitrarily. We 
choose to use a modified form of the Chebyshev polynomial of the first 
kind. In particular, we set: 


fi 


To(a) 


Pk(a) = / 


Tk(a) 




Tk(a) ^ 



k = 0 

k 2 k odd 
k > 2, k even . 


(18) 


Note that the P|^(o) conform to the required integral property; further, 
since PqCo) is identically constant with depth, it represents the depth- 
averaged component of the field. 

Further properties of the Chebyshev polynomial basis functions are 
summarized in the Appendix. 

(4b) Discretization in x and y 

Adopting a complex Fourier expansion in x, and a finite-difference 
representation in y, (17) becomes; 


(u,v,p,w) 
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0 < 1 < L 
0 £ m £ M 
0 < n < N 


(19) 


where, for instaece, = u(x,, t). The locations of the 

equivalent "grid-points" associated with this expansion are given by 
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m 


m / M 
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and 


cos [it (N-n) / N] 


It should be noted that the assumed reality of the fields (u,v,p,w) im- 
plies the associated conjugate symmetry condition: 

C“j> 'j- ‘>j- "jU * tu.j. v_*. p_*. . 


(4c) Evaluation of Spatial Derivatives and Vertical Integrals 
In representation (19), the semi-discrete form of the dynamic 
equations (10-13) can be written; 

- f» = -(us^u * vy * ^ p^[h - (l-c)»«^h 1 - ^ 

* ^ *0 
h 

= A 


+ fu = -(u«^v + vy + ws^v) - lj(p)] - (l-p)pyj - ^ 


6 (v6 v) 

» B 


( 21 ) 


if . -(U«,p ♦ vy * wy) + ^ S^(,c s^p) (22) 

= c 

and ^ 

where it is understood that the equations are to be evaluated at the grid- 
points (Xp y^, 0 ^); or alternately, in horizontal gridpoint/vertical mode 


space. 
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Here, « , 6 and { , and I^(p} denote the du^.ece forms of the dif- 

A y cr Q 

ferential and integral operators occurring in the continuum equations (10- 
13). They are evaluated as follows: 

(i) Derivatives in the x-direction are computed by direct Fourier 
synthesis using a one-dimensional discrete Fast Fourier Trans- 
form (FFT). Thus, a forward transform, followed by a multi- 
plication of the resulting Fourier spectrum by (-2nij/L ) and 
lastly, an inverse transform yields the desired x derivative. 

(ii) 6^; Derivatives in the y-direction are computed using the tradi- 

tional centered, second-order finite-difference approxima- 
tion, except at the edges of the domain (the channel walls) 
where an uncentered (one-sided) second-order formula is used. 

and (iii) and I^(p): 
y a 

Both derivatives and integrals in a can be implemented by 
constructing the appropriate matrix transformation - i.e., 
for instance, the matrix which transforms a given vertical 
column of point values into the values of the derivatives at 
those points. These transformation matrices are easily con- 
structed by reference to the known properties of the Cheby- 
shev basis functions P|^(®)* The structure of these matrices 
is discussed in the Appendix. 

(4d) Time-Stepping: Depth-Integrated Transport 
Taking k = 0, equations (20-23) describe the time evolution of the 
depth-integrated flow field. In particular - since W^jj,Q(t) = 0; i.e., 
the depth-averaged vertical velocity vanishes - the k = 0 transport field 
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horizontally non-divergent by (23). Taking advantage 
of this property, we introduce the horizontal mass transport streamfunc- 
tion where 

[u(t), v(t)]^jj^Q = «x'^^lm ’ 

Taking the curl of the momentum equations (20) and (21) in the usual way, 
we obtain a vorticity equation for the depth-integrated flow field: 

It 'Im " it Vy’’’ ' W <‘x^>*x ■ F < V* Vim = 

*'lni C*x 1'y - V *x * 'x®0 ’ Vo^lm 

Note that (Aq)^i^ and (8 q)ijj| are the depth-integrated components of the 
cumulative right-hand side terms in equations (20) and (21). This compoiT- 
ent is obtained from the associated physical space (vertical) distribu- 
tion of values by a simple matrix transfonnation - see the Appendix. Not 
only does the introduction of ijj automatically guarantee horizontal non- 
divergence, but it eliminates the pressure gradient term associated with 
the (unknown) reaction pressure supported by the rigid lid. A solution 
to equation (24) is fully prescribed by specifying the (constant) bound- 
ary values of on the channel walls - that is, by setting: 

<jilo(t) = 'i'Q(t) and ,|,^,^(t) = f|^(t) . (25) 

The vorticity equation (24) is advanced in time using a second-order 
Adams-Bashforth approximation. The resulting fully discrete equation is: 
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where q is the time step index. By applying a forward transform in x, the 
i'elmholtz-1 ike equation (25) ~ and associated boundary conditions (25) - 
can be reduced to a sequence of (M-2) x (M-2) complex tri agonal matrix 
equations, one for each Fourier component (j) in the x expansion. These 
triagonal equations can be efficiently solved to yield the updated stream- 
function field 

(4e) Time-Stepping: Internal Velocity Modes (1 £ k £ N-2) 

The (N-2) lowest order internal modes of the velocity distribution - 
that is, for 1 £ k £ N-2 - can be obtained by direct time-step- 

ping of equations (20 and (21). Using an Adams-Bashforth time-step on 
the right-hand side terms A and B, and a Crank-Nicolson approximation on 
the Coriolis terms, we get: 
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where the indices 0£l£L, 0£m£M and 1 £ k £ N-2. The straightfor- 
ward algebraic solution to these equations is: 
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Hmk 
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where 


o = fAt/2 . 


{4f ) Time-Stepping: Density Equation (0 < k < N-2) 

The 0 £ k £ N-2 vertical modes of the density field are advanced in 
time using equation (22), Without further elaboration, the resulting 
fully discrete form of (22) is: 


q+1 q 
‘’imk “ ‘'imk 


+ (4i-) 

' t ' ^Imk 


- (|^) C 


g-1 

Imk 


(31) 


where 0£l£l, 0£m£Mand0£k£ N-2. 

(4g) Application of Boundary Conditions Using the Tau Method 
Note, in Sections (4d, e, f) that only the lowest (N-1) modes of 
(u,v,p)*^^^ have been obtained by time integration of the dynamic equa- 
tions (20-22). The remaining two highest vertical modes - that is to say, 
(u,v,p)^^^ with 0£l£l, 0£m£M and N-1 £ k £ N - are chosen such 
that the resulting (u,v,p)*^^^ fields identically satisfy the prescribed 
boundary conditions on a = £ 1 (equations 15-16). This is the so-called 
tau method, and can be shown to deliver solutions of comparable accuracy 
to those obtained using Galerkin or collocation approximation techniques. 
Definitions, further discussion, and sample comparisons of these spectral 
approximation techniques has been given by Gottlieb and Orszag (1977). 

Having obtained for the lower vertical modes 0£ k £ N-2, 

the solution of a simple algebraic problem at each horizontal location 
(Xii^, y^j^) provides the required values of (o,v,p)^iJjJ for k = N-1, N. The 
precise structure of these algebraic equations is discussed in the Appendix. 
(4h) Determination of the Vertical Velocity Field 
Having obtained a complete specification of the u,v, and p fields at 
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time step (q+1) by the methods outlined above, the vertical velocity field 
at the interior vertical gridpoints (1 < n <_ N-1) can be computed: 

'^?mn irr ^ ^y ^*^''lmn^^ 

lin n 

where 0^1 <L, 0£m<M and 1 £ n < N-1. In the sigma coordinate sys- 
tem, = 0. 

5. Model Performance and Evaluation 

The numerical model described in the preceding sections has been 
validated on a variety of simple test problems for which analytic 
solutions are available. Results for one such representative class of 
tests are given next. 

With uniform rotation and density (f = fg and p = 1), and with the ex- 
ponential bottom profile h(x,y) = Hge“^®-y, nonlinear topographic wave 
solutions of equations (20-24) can be found. The solutions are depth- 
invariant; hence, they can be written explicitly in terms of the trans- 
port streamfunction. The analytic solution is; 

♦a(x,y,t) = e-“> sin (.jy/L^) , (33) 

with the accompanying dispersion relation 

0) = - 2fgak / (k^ + + a^). (34) 

The transport components (hu, hv) can, of course, be recovered from the 

solution for by the diagnostic relations (uh, vh) = (-ip , 4 > ). 

y X 

For the topographic wave test problems, the following nominal values 
for the dimensional parameters have been used; 
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(i) 

along-channel length 

B 10® cm B 1000 km 



(ii) 

cross-channel length 

■ 5 X 10^ cm = 500 km 



(iii) 

depth at y = 0 

Hq = 5 X 10^ cm = 500 m 



(iv) 

Coriolis parameter 

fQ = 1.0 X 10”^ sec“^ 

I:! 

and 

(V) 

constant density 

p = 1 gm cm~® 


The choice of a specific analytic solution \p^ (x,y,t) involves the speci- 
fication of three physical parameters . These are: the along-channel wave- 
number, k; the cross-channel mode number, j; and the topographic e-folding 
scale length, o. Having chosen these parameters, the frequency of the wave 
follows from (34). 

In ascertaining the accuracy with which the numerical model recovers 
the nonlinear wave solutions (33), the important non-dimensional computa- 
tional parameters are: 

(i) (Number of points in x/x wavescale) = L/k 

(ii) (Number of points in y/y mode scale) = Mk/j or M/oL^ 
and (iii) (at / wave period) = atw/2ir . 

Since the analytic solution is depth-invariant, there is no obvious 
meaningful measure of the formal accuracy of the a expansion. However, it 
will be important to note any spurious generation of higher modes (k > 0) 
by the numerical model; since these modes identically vanish for the 
analytic solution, we hope that they remain small in the model solutio . 

Table I shows the results of these analytic tests. As well as the 
relevant parameter values, the table lists the normalized RMS differences 
between the model -generated and analytic wave solutions for \ii, u and v 
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after an elapsed time of one wave period. Here, the normalized RMS error 
in t|) is defined as 


RMS (i|^) 



n 2 / 

<♦ ■ *4)1™ / 

msi 



M 

E 

m^l 




with similar definitions for RMS(u) and RMS(v). 

For these wave problems and the range of parameters studied, RMS er- 
rors are generally about one percent after one wave period. Model errors 
depend most sensitively on the* non-dimensional time-step. (Compare, for 
instance, cases 4 and 5. No<.<i the nearly quadratic dependence of RMS er- 
ror on at.) The relative insensitivity of RMS error to the effective x 
resolution is a reflection of the high accuracy of the Fourier expansion 
technique. In particular, it is well known that Fourier simulations of 
wave propagation are free from any phase error attributable to space dif- 
ferencing for any resolved wave (Orszag and Israeli, 1974). While this is 
certainly not true of the second-order finite-difference representation 
in the y-direction, the analytic solutions examined here have phase prop- 
agation in X only. Since the cross-channel modal structure is invariant 
in time, RMS errors are only weakly dependent on the cross-channel reso- 
lution parameters. An examination of cases 1, 4 and 6 indicates, however, 
that ay-related RMS errors comparable to those arising from time-differ- 
encing begin to appear for these wave solutions when (M/j) approaches 
0(10). (Compare also cases 7 and 8.) Lastly, we note that the simulated 
wave solutions remain highly depth-invariant througnout the experiments. 
Higher mode amplitudes are typically smaller by four orders of magnitude 
(nearly machine accuracy) than the barotropic mode. 
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These topographic wave tests have been run at the Woods Hole Oceano- 
graphic Institution on a VAXll-780 system. Calculations were carried out 
primarily in-core, with a minimum of "paging". The required cpu time per 
mroel time step is shown in Table II as a function of spatial resolution. 
These figures indicate that execution time increases roughly linearly with 
the number of horizontal grid points, and perhaps slightly more rapidly 
with the number of vertical expansion functions. To put these cpu times in 
an appropriate perspective, a two-month simulation using a (32 x 49 x 5) 
grid, and a time step of 0.1 day would require roughly 3,5 cpu hours on 
the VAX, and much less on a more powerful machine. Even on a machine like 
the VAX, however, this cpu requirement is modest. 

Prototype test problems characterized by non-trivial vertical struc- 
ture have also been run, with roughly comparable error levels, "^cken to- 
gether, these results indicate that the four-dimensional model described 
here can be used to provide accurate and efficient numerical exploration 
of coastal/deep ocean coupling and physics. 
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TABLE I 

Normalized RMS Model Errors in Streamfunction and Velocity (u,v) 
After One Period for Nonlinear Topographic Wave Solutions 


Case 

L^o(M/(iL^) j(M/j) 

k(L/k) 

L 

M 

N 

atw 

RMS«I- 

(*/.) 

RMS u 
(•/•) 

RMS V 
(•/.) 

1 

1.0 (49) 

1 (49) 

4 (4) 

16 

49 

5 

0.01 

1.07 

1.06 

1.07 

2 

1.0 (25) 

1 (25) 

4 (4) 

16 

25 

5 

0.01 

1.14 

1.14 

1.14 

3 

1.0 (49) 

1 (49) 

4 (8) 

32 

49 

5 

0.01 

1.19 

1.33 

1.61 

4 

1.0 (49) 

2 (25) 

4 (4) 

16 

49 

5 

0.01 

1.26 

1.26 

1.26 

5 

1.0 (49) 

2 (25) 

4 (4) 

16 

49 

5 

0.02 

4.46 

4.41 

4.46 

6 

1.0 (49) 

4 (12) 

4 (4) 

16 

49 

5 

0.01 

2.94 

3.08 

2.94 

7 

2.0 (25) 

2 (25) 

4 (4) 

16 

49 

5 

0.01 

1.34 

1.33 

1.33 

3 

2.0 (13) 

2 (13) 

4 (4 

16 

25 

5 

0.01 

2.19 

2.53 

2.17 

9 

1.0 (25) 

2 (13) 

4 (4) 

16 

25 

10 

0.01 

1.88 

2.30 

1.88 





TABLE II 

Required cpu Time Per Model Time-Step on a VAXll-780 
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APPENDIX 

DETAILS OF THE CHEBYSHEV-TAU SOLUTION TECHNIQUE 


Recall from Section 4 that the vertical (a) dependence of the depen- 
dent variables (u,v,p,w) has been represented as an expansion in the mod- 


ified Chebyshev basis function P|^(o), where 




To(») 

o 

II 


TflCa) 

k > 0, odd 

(A.1) 

,3 

k > 0, even 



- 1 


and T|^(o) is the k— degree Chebyshev polynomial of the first kind. Typi- 
cally defined over the interval -1 < o <. Tj^(o) = cos[k cos~^(o)] where 
-ir _< cos“^(o) < 0. For the purposes of the following discussion, we in- 
troduce an arbitrary one-dimensional variable n, and set 


= n 


N 

(a ) = Z 

" k=0 


\ 'k 




0 < n < N 


(A.2) 


The "collocation" points or equivalent gridpoints, o^, are conveniently 
chosen to correspond to the extrema of T|^(a); hence, 

c = cos [ir (n-N)/N] . 0<n<N (A.3) 

Given n(a^) at the (N+1) collocation points (0 < n N), the equiv- 
alent expansion coefficients ri|( (0 £ k <_ N) can be recovered in several 
ways. For instance, a Fast Cosine Transform (FCT) can be used - see, 
e.g., Haidvogel and Zang (1979). Alternately, a linear (matrix) trans- 
formation can be constructed to compute the n|^ from Pp. The required 
transformation matrix F is given by 
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F 


Pq(-i) Pi(-i) .... P^C-i) 
Pq(®i) Pi(°x^ *••• 

yi) Pi(i) .... p„(i) 


(A.4) 


Considering n„ and n|< to be column vectors, it is easily seen that 



'""k * "n 

0 < n, k < N 

(A.5) 

and, hence, that 

C-1 

^ '’n - \ * 


(A.6) 


The identical approach can be used to form matrix operators which 
perform a vertical derivative or integral. For instance, 

^‘°n* = ' 0 % = “"n (A.7) 


where 


0 


’po(-i) y-1) .... Pn(-i) 
y°i) picoj) .... Pn(oi) 



I Pfld) 



(A.8) 


and 


’’k*") ' '’kt.) = ^ T|^(o) 


(A.9) 


n sin (ne) 
sin & 


0 = cos“^(o) 



Likewise, 


(A.IO) 
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(Aai) 


and 



I (1 - o^) k . 0 

I (1 - 4 ) k . 1 


cos (n+l)x 

COS (n-l)x ' 

0 

1 1 

2 (n+1) 

■ "i (n-1) 

< 

1 cos"^(o ) 
0 " 

/cos (n+l)x 

[~rmr 

COS (n-l)x ^ 

- 

1 i 1 

cos {a) 
n 

* (1 - On) 

1 


(k^ -1) 



(A.12) 


k > 1, odd 


k > 1, even 


As described in Section 4, the tau approximation replaces the dynamic 

equations for k = N-1, N - that is, those which dictate the dynamical 

evolution of the highest two vertical modes - with the two boundary 

conditions at top and bottom (o = + 1). The implementation of the tau 

method proceeds as follows. Suppose that our time-stepping scheme has 

provided values of (0 1 1 N-2), but that our boundary conditions 

require that 


n (-1) = 0 
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and 

If(i) 

» Y 

Then, 

Vii-i) v-i)' 


"nq+l" 

Vi 


Vi(i) pjd) 


“n 


where 


\ = 




k = 0, .... N-2 
k = N-1. N 


- n{-l) 1 


dfi 

Tr 

da 


U) 


(A.13) 


(A. 14) 


The algebraic equation (A.13) is readily solved. Note that the same seq- 
uence of steps leading to (A.13, A.14) yield similar simple expressions 
for (k = N-1, N) for various alternate top and bottom boundary 


conditions. 
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